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We propose a comprehensive framework for flexible additive re- 
gression models for correlated functional responses, allowing for mul- 
tiple partially nested or crossed functional random effects with flexi- 
ble correlation structures for, e.g., spatial, temporal, or longitudinal 
functional data. Additionally, our framework includes linear effects of 
functional covariates and linear or smooth effects of scalar covariates 
that may vary smoothly over the index of the functional response. It 
accommodates densely or sparsely observed functional responses and 
predictors which may be observed with additional error. Inference in 
this framework can be performed by standard software for generalized 
additive models (GAMs), allowing us to take advantage of established 
robust and flexible estimation algorithms. We provide easy-to-use 
open source software for the specific purpose of fitting functional ad- 
ditive mixed models in the pf f r () function for the R-package refund. 
Simulation experiments show that the proposed method recovers rel- 
evant effects reliably, handles small group sizes and/or low numbers 
of replications well and also scales to larger data sets. Application ex- 
amples on spatial and longitudinal functional data demonstrate that 
the proposal yields flexible model specifications that do justice to 
complex data situations and yield interpretable results. 



1. Introduction. In recent years, many scientific studies have collected 
functional data that exhibit correlation structures amenable to explicit mod- 
eling. Such structures may arise from a longitudinal study design (e.g. Morris 
and Carroll, 2006; Greven et al., 2010; Zipunnikov et al., 2011; Goldsmith 
et al., 2012) or spatial sampling of curves (e.g Giraldo et al., 2010; Delicado 
et al., 2010; Nerini et al., 2010; Staicu et al., 2010). Simultaneously, regres- 
sion for independent functional responses (Faraway, 1997) has made large 
advances, including both multiple scalar (Reiss et al., 2010) and multiple 
functional predictors in concurrent or more general relationships (Ivanescu 
et al., 2011). 

Our work is motivated by a longitudinal neuroimaging study on multiple 
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sclerosis (MS). MS is an autoimmune disease of the central nervous system 
(CNS) that can lead to severe patient disability. It damages white matter 
tracts in the brain due to lesions, axonal damage and demyelination. Diffu- 
sion tensor imaging (DTI) is a magnetic resonance imaging (MRI) technique 
that is able to resolve individual white matter tracts in the central nervous 
system (Basser et al., 2000; Gossl et al., 2002), and is thus a very useful tool 
in monitoring disease progression in MS patients. Our data set contains re- 
peated measurements for 106 MS patients and 42 healthy controls. At each 
visit, DTI was carried out and fractional anisotropy (FA) was determined 
along some of the most important white matter tracts in the brain. FA de- 
scribes the relative degree of anisotropy of diffusivity along white matter 
tracts. It may be decreased in MS patients and thus serves as a marker of 
disease progression here. FA values were extracted along the corpus callo- 
sum (CCA, connecting the left and right hemispheres of the brain) and the 
corticospinal tract (CST, connecting the brain and the spinal cord), yielding 
two curve measurements per visit and subject. The goal of our analysis is 
to evaluate the prognostic power of tract profiles at a previous visit for the 
next visit in MS patients, while accounting for effects of age, gender and 
the elapsed time between visits as well as within-subject correlation of the 
repeated measurements. 

There is a need for regression models for functional responses that can 
accommodate more general correlation structures suitable for longitudinal or 
spatial data via functional and scalar random effects. At the same time there 
is a need for a more general modeling framework for functional responses 
that accommodates flexible effects of both scalar and functional covariates. 
In this paper, we generalize the methods described in Ivanescu et al. (2011) 
to allow for correlated responses. Our approach is able to handle densely or 
sparsely observed functional responses and predictors which may be observed 
with additional error. We introduce additive mixed models for functional 
outcomes that inlcude both functional and scalar covariates and describe 
the estimation and inference machinery. Specifically, we consider models of 
the general form 

Yij(t) = a(t) + / Xij(s)P(s,t)ds + ji(zuj,t) + 72(221.7) + z 3i jS 3 + 
(1) Js 

zujSi(t) + b 0i + buuuj + B Qi (t) + Bufyuuj + eij(t), 

with functional responses observed over a domain T for subjects or 

groups i = 1 , . . . , M and repeated observations j = 1 , . . . , rij . A total of 
n = n i functional observations are available. The additive predictor can 
contain any combination of 
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- a smooth functional intercept a (t), 

- linear effects /3/(s, t), f = 1, . . . , F of functional covariates Xf(s), with do- 
mains Sf that may or may not be the same as the domain of the functional 
response, 

- functional or constant smooth effects of scalar covariates z r , r = 1, . . . ,p: 

j(z r ,t), j(z r ), 

- functional or constant linear effects of scalar covariates z r , r = l,...,p: 

- scalar random effects b q i with associated scalar covariates u q , 

- functional random effects B q i{t) with associated scalar covariates u q . 

Note that smooth terms like 7(2;, t) or 7(2:) also accommodate vector-valued 
covariates z like latitude-longitude pairs for spatial coordinates. For simplic- 
ity of presentation, we include only one term of each kind, but the extension 
to several terms is obvious and available in our implementation. 

We assume that a(t), (3f(s,t), f = 1,...,F, j r (z, t), j r (z) and 8 r (t), 
r = l,...,p are smooth but unknown functions, possibly dependent on 
several covariates. Error terms tij{t) are assumed to be independent and 
identically distributed (i.i.d.) Gaussian variables with mean zero and con- 
stant, unknown variance across T, such that e(t) can be viewed as a white 
noise error process, bj = (boi,bu) are assumed to be i.i.d. mean-zero Gaus- 
sian random vectors with general unknown covariance matrix D. Functional 
random effects B q i(t) are modeled as realizations of a mean-zero Gaus- 
sian random process on T x { 1 , . . . , M} with a general covariance function 
K Bq (t,t' = Cov (B q i(t), B q i/(t')). Different functional random effects 
B q i{t),B q 'i(t) are assumed mutually independent. Note that our model class 
and software also admit multiple partially or completely nested or crossed 
grouping factors for both scalar and functional random effects. The smooth- 
ness assumptions on all model components besides €ij(t) ensure smoothness 
of Yij (t) up to the white noise measurement error €ij (t) . 

Estimation and inference in our approach are based on a representation 
of model (1) in terms of a conventional scalar-on-scalar penalized additive 
regression model. The main advantage of this strategy is that we gain imme- 
diate access to the powerful and versatile inference machinery developed for 
this broad class of models over the last years. A second major advantage of 
large practical importance lies in the fact that this representation naturally 
accommodates noisy functional responses Y(t) observed on irregular and/or 
sparse grids. 
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8s To the best of our knowledge, our proposal is the first description and im- 

86 plementation that allows such a high level of flexibility for a regression model 

87 for functional responses: Previous approaches either limit the predictor to 

88 the effect of a single functional covariate and a functional intercept, such as 

89 the linmod function in package fda (Ramsay et al., 2011) for the R language 

90 (R Development Core Team, 2011) or to scalar covariates only, such as the 

91 fosr function for scalar-on-function regression (Reiss et al., 2010) in the R- 

92 package refund (Crainiceanu et al., 2011). Like a simpler precursor of our 

93 approach described in Ivanescu et al. (2011), neither of those approaches is 

94 able to account for correlated functional responses. Our proposal has some 

95 similarities with the regression models developed for independent or longitu- 

96 dinal scalar responses in Goldsmith et al. (2011) and Goldsmith et al. (2012), 

97 since we also base our inference on additive mixed models for scalar-on-scalar 

98 regression and take the same approach to deal with noisy, potentially sparse 

99 functional predictors. However, the extension to functional responses and 

100 functional random effects with flexible correlation structure presented in 

101 the following is non-trivial. 

102 Functional principal component approaches for longitudinal (Greven et al., 

103 2010), multilevel (Di et al., 2009) or spatially correlated (Staicu et al., 2010) 

104 functional data are primarily concerned with decomposing observed variabil- 
es ity of detrended functional data Y%j{t) — E(Yij(t)), but are not aimed at flex- 

106 ible modeling of the mean structure. Previous spline- or wavelet-based func- 

107 tional linear mixed models such as Guo (2002); Morris and Carroll (2006); 

108 Antoniadis and Sapatinas (2007) are limited to functional linear effects of 

109 scalar covariates and responses observed on a regular grid and are unable 
no to deal with (non-concurrent effects of) functional covariates. Going beyond 
in the model of Brumback and Rice (1998), our approach accommodates func- 

112 tional random effects with between-subject correlation structure and offers 

113 much more flexibility for the model's fixed effects. 

114 The paper is organized as follows: Sections 2 and 3 develop our general 
us approach for penalized regression for correlated functional responses and its 

116 estimation in a mixed model framework. Section 4 discusses identifiability 

117 issues. Our method is evaluated in a simulation study, in an application to 
us the motivating longitudinal DTI study and the spatial Canadian Weather 

119 data set in Section 5. Section 6 closes with a discussion and outlook. 

120 2. Penalized regression for correlated functional data. As noted 

121 by Ramsay and Silverman (2005, p. 259), models for functional responses 
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122 are very similar to varying coefficient models. We exploit this similarity in 

123 order to re- formulate estimation and inference for model (1) in terms of 

124 estimation and inference for a conventional additive mixed model (AMM) 

125 in which effects vary smoothly over the responses' domain T ■ This allows us 

126 to fit complex and versatile models for functional outcomes and covariates in 

127 realistic data situations using the flexible and highly-developed estimation 

128 and inference routines available for a very broad class of AMMs in the mgcv 

129 package (Wood, 2011) for R. In the following, we drop the double index for 

130 observations and groups to simplify notation where possible. For a given 

131 matrix Z, we denote the i-th row by Z%. and the element in row i, column 

132 j by Zij. 



2.1. Representing regression models for functional responses as varying 
coefficient models. Functional responses Yij(t) are observed on a grid of Gij 
evaluation indices = (Uji, ■ ■ ■ >%Gj,) ■ hi order to simplify notation and 
without loss of generality, we assume identical grids tij = t = (ti, . . . ,tc) T 
Vi,j in the following, but note that functional responses Y(t) observed on 
irregular and/or sparse grids are naturally accommodated in the rephrased 
model formulation given in (2) below. We write 



Y n (h) 



Y lni (h) 



YMn M {h) 



Y n (t G ) 



Yi ni (t G ) 



Ymu m (*gO 



vec(Y 



133 for the n x G matrix of observed response values and y 

134 (yJi, • • • , yj ni , • • • , VM nM ) T f° r tli e vector that holds the Gn concatenated 

135 function evaluations yij = (yij 1 , ■ ■ . ,Uij G ) T which form the rows of Y. In 

136 doing this, we switch from a functional data perspective, which views each 

137 Yij(t) evaluated on £ as a single observation, to a longitudinal data per- 

138 spective, which views the function values observed at timepoints £ as a vec- 

139 tor of separate responses (yij 1 , ■ ■ ■ ,Uij G ) T ■ The smoothness assumption on 

140 E(Yy(t)) is preserved implicitly by enforcing smoothness across T for a(t) 

141 and the various (3f(s,t), S r (t), j r (z r ,t) and B q i(t). 

142 In this notation, model (1) can be equivalently expressed as: 



Vijg = «(*s 



Xij(s)f3(s,tg)dS + Jl(zUj,tg) + j2(Z2ij) + 



(2) JS 

Z3lj^3 + Z±ij8±{tg) + + hiUlij + B 0i (tg) + Bli(tg)UUj + Cijg. 
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143 The assumption of white noise error terms translates to e^g '~ ' iV(0, a 2 ). 

144 In this perspective, linear functional effects zS(t) of scalar covariates are 

145 simply smoothly varying coefficients of t and smooth functional effects j(z, t) 

146 are bi- or multivariate functions of a scalar or possibly vector- valued covari- 

147 ate z and the index of the response. The extended covariate matrix for model 

148 (2) is constructed by repeating each entry in the original covariate vectors 

149 G times as described in Ivanescu et al. (2011). We also add a vector of n 

150 repeats of t as the response index variable t. In the remainder of this paper, 

151 we assume that the appropriate concatenation and repetition operations de- 

152 scribed above have been performed and let our notation refer to the extended 

153 dataset. 

154 If the functional covariates Xf(s) are observed on incomplete or sparse 

155 grids a preliminary step that allows reconstruction of the curves on a regular 

156 and dense grid is required. Techniques such as functional principal compo- 

157 nents analysis (FPCA) (Yao et al., 2005), longitudinal FPCA (Greven et al., 

158 2010), or spline smoothing (James, 2002) can be used. This is similar to the 

159 case of scalar-on-function regression developed in Goldsmith et al. (2011, 
wo 2012). 

161 In the following sections, we describe the modeling methodology for the 

162 various effects of the scalar and functional covariates of model (1). Our ap- 

163 proach is based on using penalized splines to model the smooth effects. Be- 
rn cause the smooth effects exhibit different degrees of complexity, we divide our 

165 presentation into three parts: modeling the smooth covariate effects (Section 

166 2.2), the integrals of smooth effects (Section 2.3), and modeling functional 

167 random effects (Section 2.4). We will drop the double index for units i and 

168 observations j to reduce notation in the following and use a single index i 

169 for all observations instead. 



2.2. Spline basis representations of smooth effects. Smooth effects of a 
single covariate (i.e., a(t), j(z), S(t)) are represented as weighted sums of 
spline basis functions, with appropriate roughness penalties on the spline 
coefficient vectors. For example, 

K t 

k t =l 

170 for i = 1, . . . , N, r = 1, . . . ,p, with a penalty term X r 6 rT P r r added to the 

171 log-likelihood. &k(zi) is the value of the /c th basis function evaluated at z„, 
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172 is the i row of the N x Kt matrix 3> r of basis function evaluations, 

173 r = (9 r ,i, • • • , @r,K t ) is the- vector of spline coefficients. P r is a known and 

174 fixed Kt x Kt positive (semi-)definite penalty matrix and A r is a positive 

175 smoothing parameter controlling the trade-off between goodness of fit and 

176 the smoothness of 7(2™ ) (Wood, 2006a, ch. 4.1). 

Smooth effects of two or more covariates (e.g. 7(2, t)) are represented as 
tensor product splines, i.e., 

K z K t 

177 for i = 1,...,N, r = l,...,p. They are associated with penalties 

178 rT (X r (P r <8> Ixt) + ^rt(lK x <S> -P 4 )) # r , where is the kxk identity matrix 

179 and P r ,P t are penalty matrices for the marginal smooth terms. As before, 

180 <&* is the g-th row of <& with entries $>k(t g ), k = 1, . . . The matrices 

181 $ 2:r and <& contain evaluations of the respective marginal basis functions, 

182 and the design matrix for the tensor product spline is given by the row- wise 

183 tensor product of their rows <&f r and (Wood, 2006a, ch. 4.1.8). This 

184 flexible construction is valid for any combination of bases with quadratic 

185 penalties. 

186 Varying coefficient terms z r S r (t g ) are correspondingly represented as 

187 z r ij& g 9 r , leading to a design matrix $ r for r that is the Kronecker product 

188 of a vector containing the z r ij and a matrix with rows 

189 Each scalar covariate z r can only be associated with either a linear vary- 

190 ing coefficient term z r 5 r (t) or a constant smooth term 7 r (z r ) or a smooth 

191 functional term 7 r (z r , t). 



2.3. Spline basis representations of integrals of smooth effects. As for the 
other smooth effects, we represent each effect surface /3/(s, t), f = 1, . . . , F 
for the linear effect of a functional covariate Xf(s) as a tensor product spline. 
The effect evaluated at s = (si, . . . , sh) T and i 9 is /3/(s, t s ) ~ <&**0^, where 
<I>g* is the -ff x K s Kt matrix of tensor product basis functions in s and t 
evaluated for all index values s and specific index value t g . The integrated 
effect can be approximated as 

H 

/ X fi ( 8 )p(s, t g )ds « 2 w h X fi (a h )P(8 h ,t g ) = (L ■ X fi ,)*fOf = 
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where • denotes elementwise multiplication and L is a 1 x H matrix con- 
taining the quadrature weights Wh for a numerical integration scheme. The 
rows of the new basis 

H, = ( L - x f*-)*f 



192 are 1 x K s Kt vectors containing the approximate integrals of the tensor 

193 product basis functions scaled by Xfi(s). Thus, the integrated effect can be 

194 represented in terms of a conventional tensor product spline with modified 

195 basis functions (c.f. Wood, 2011). Its penalty is given by the standard tensor 
we product penalty 0f T [Xf s (P s ® Ix t ) + \ft(Iic s ® P*)) & ■, where I k is the 

197 k x k identity matrix and P s , P l are penalty matrices for the marginal 

198 smooth terms in s and t. 



2.4. Spline basis representations of functional random effects. Functional 
random effects B q i(t), q = 0, . . . , Q are represented as smooth functions for 
each observational unit i = 1, . . . , N and modelled using penalized splines 
with a tensor product design matrix and a tensor product penalty. For sim- 
plicity of notation, we assume a single grouping variable in the following, 
but extensions to multiple (partially) nested or crossed random effects are 
trivial. The design matrix for functional random effects is given by the Kro- 
necker product of the nxM incidence matrix R for the levels of the grouping 
variable and a G x Kt spline basis matrix <&* with Kt basis functions eval- 
uated on t, leading to a N x MKt design with Kt basis functions for each 
grouping level i = 1, . . . , M: 

M K t 

m=l k=l 
M K t 

B u (t g )uij ~ Y ^2s im u ij B k (t g )e mk = mj ■ (Ri. (8) <S> g .)0 l = &\ g l , 

m=l k=l 

199 with Kronecker's delta 5i m . 

200 Smooth random effect functions B q i(t), B^i (t) for different units i, i 1 share 

201 common smoothing parameters. The penalty matrix 

202 \{P q ® IkJ + Aqt(JM ® P l ) is a weighted sum of penalty matrices, where 

203 P q is a M x M precision matrix for the dependence structure between units 
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204 that encourages similarity between the units' spline coefficient vectors and 

205 .P* is a Kt x Kt smoothness penalty matrix that is applied to every subvec- 

206 tor of unit-specific spline coefficients. The smoothing parameter \ q controls 

207 the relative precision of the inter-unit variability, i.e., it controls the degree 

208 of similarity between the estimated functional random effects for the dif- 

209 ferent group levels, while X q t is a global smoothing parameter controlling 

210 the common roughness of the functional random effects. To see this, con- 

211 sider that the penalty for the subvector (Oik, ^2jfc, • • • > ^Mfc) of coefficients 

212 for the k-ih basis function <£&(£) for group levels i = 1, . . . , M is given by 

213 X q P q + \ q tPjl k , i.e., the precision matrix specifiying inter-unit variability 

214 plus the precision for the k-ih basis function, and that the penalty for the 

215 subvector (0 m i, 0m2, ■ ■ ■ , QmK t ) of coefficients for the spline basis functions 

216 in the m-th group level is given by \ q Pmm + Agt-P*, i-e., the precision for the 

217 m-th group level plus the smoothness penalty matrix in i-direction. 

218 If observational units on different levels of the grouping factor are assumed 

219 independent, P q is simply the identity matrix with as many columns/rows 

220 as there are levels of the grouping factor, i.e. P q = Im- Much more gener- 

221 ally, P q can be used to incorporate a dependence structure between units: It 

222 can be a (partially improper) precision matrix of a random field with known 

223 correlation structure, implied, for example, by the spatial or temporal ar- 

224 rangement of the units, such as a Gaussian Markov random field (GMRF) on 

225 geographical regions for conditionally auto-regressive (CAR) model terms. 

226 Alternatively, (P q )~ 1 can be defined using any valid correlation function 

227 based on spatial, temporal, or genetic distances between levels of the group- 

228 ing variable, for example. If the grouping variable is simply an index of the 

229 observations, this construction can also be used to model correlated func- 

230 tional residual curves to account for spatial or temporal autocorrelation of 

231 the responses. This innovative definition of functional random effects ad- 

232 mits very flexible model specification, since any combination of spline basis, 

233 smoothness penalty and between-subject correlation can be used for func- 

234 tional random effects, allowing, for example, for spatially correlated func- 

235 tional residuals with periodicity constraints for the Canadian Weather data 

236 (see Section 5.2 and Appendix C.2 of the online supplement). 

3. Mixed model representation of a GAM. Consider model (1) 
and recall the spline basis representations of the various smooth effects given 
in the previous Section 2. Collecting all the pieces, model (2) can be re- 
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written as 



(3) yijg = &ij 9 + e ijg ; e ijg ~ N(0, a. 



e )i 



where <f> — \<f>®- ■ ■ ■ ■ <b p ■ &> l ■ . . . ■ <6^ ■ &9. • . . . • d>9. 1 

237 is a 1 x d-row vector concatenating the ijg-th. rows of 

238 - the design matrix 3>° for the functional intercept a(t) (Section 2.2), 

239 - (tensor product) design matrices <& r , r = 1, . . . ,p for the various smooth 

240 or linear effects of scalar covariates (Section 2.2), 

241 - integrated scaled B-spline tensor products for the linear effects of func- 

242 tional covariates Xf(s), f = 1, . . . , F (Section 2.3), 

243 - and the design matrices q = 0, . . . , Q for the functional random effects 

244 (Section 2.4). 

245 The d x 1 vector 6 contains all the stacked coefficients associated with the 

246 respective design matrices. 

247 Since spline basis dimensions have to be sufficiently large to ensure enough 

248 flexibility, maximum likelihood estimation of model (3) will lead to substan- 

249 tial overfitting. This can be avoided by using a penalized likelihood approach, 

250 where the penalties are designed to suppress excessive wiggliness of the var- 

251 ious components. We use quadratic penalties for univariate smooths, and 

252 tensor product penalties for bivariate smooths (Wood, 2006a,b), as described 

253 in Section 2. For notational simplicity, we re-denote the various penalty ma- 

254 trices and smoothness parameters following the notation of Wood (2006a) 

255 by assigning a sequential index i = 1, . . . , 1 + p + F + (Q + 1) for all the 

256 smooth terms in model (3). Let Qi,\i and Si denote the coefficient vector, 

257 smoothing parameter and penalty matrix respectively corresponding to the 

258 I th univariate smooths, and let 0£, A^i, \£2 and $n, &£2 denote the coefficient 

259 vector, the two smoothing parameters and tensor product penalty matrices, 

260 respectively corresponding to the £ th bi-variate smooths. Next, we re- write 

261 the roughness penalties above in terms of the full coefficient vector so 

262 that a general quadratic penalty 6j S(0( becomes 6 T S^O, where Sg is just 

263 Si padded with zero such that 9 T SiO = 6j SpO?. Likewise a general ten- 

264 sor product penalty Oj(\aSii + A^S^)^ becomes \nO T S^O + Xf20 T St3O 

265 where Sn and Si2 are just Sa and Si2 respectively padded with zeros such 

266 that eJ(X ei Sii + \nSn)0l = XaO T SaO + X^S^O- Therefore the penal- 

267 ized likelihood criterion to be minimized is reduced to: 



a~ * — ' ' — ' c; 

1,0 1=1 
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268 where I is used to index all the smoothing parameters A; and the correspond- 

269 ing penalty matrices Si and is the G x d matrix obtained by column 

270 stacking the Bij g over g = 1, . . . , G. The total number of smoothing param- 

271 eters is L = 1 +p\ + 2p2 + 2F + 2{Q + 1), p = p\ + P2, where p\ denotes the 

272 number of univarate smooths of scalar covariates (i.e. 5 r (t)z r , J r (z r )) and p2 

273 denotes the number of bivarate smooths of scalar covariates (i.e. "f r (z r ,t)). 

274 The F linear effects of functional covariates and the Q + l functional random 

275 effects have two smoothing parameters each. 

2 

276 Let ti = ^ and use similar arguments as in Wood (2006a), to obtain the 

277 solution of (4) as the best linear unbiased predictor in the mixed effect 

278 model (MEM): 



279 where S~ denotes the generalized inverse of S, and N(0, S~) is an improper 

280 Gaussian prior with positive semi-definite covariance matrix S. The prior 

281 impropriety results from rank deficiency of X^z 7 "; -1 ^ anc ^ can ^ e overcome 

282 by considering the eigen-decomposition of ^2iT^ 1 Si = UDU T where U is 

283 an orthogonal matrix of eigenvectors, D is a diagonal matrix of eigenvalues 

284 with c zero-eigenvalues and the partition U = [Upo ■ Urq], where Ufo has c 

285 columns. These are well known issues in the literature on penalized regression 

286 splines (see Wood (2006a,b)) and we discuss them briefly for completeness. 

287 Therefore if Vij = <&ijUpo, Zij = $^17ro and Si = U^SiIIrq then model 

288 ( 5) can be re-written using a more convenient MEM representation as 



290 formed coefficients which are penalized from those which are not; the prior 

291 for Oro is now a proper prior. 

292 One of the main advantages of formulating the penalized likelihood op- 

293 timization as estimation in an MEM is that the smoothing parameters A; 

2 

294 can be treated as variance component parameters (A/ = ^f) and thus can 



(5) 






289 where (^Q'^ro) T ^ s a re-parameterization of which separates the trans- 
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295 be estimated using restricted maximum likelihood (REML). In particular, 

296 Reiss and Ogden (2009) have shown that smoothing parameter selection 

297 with REML outperforms other selection criteria such as Akaike's informa- 

298 tion criterion or generalized cross-validation. Additionally, the inferential 

299 machinery developed for mixed models can be used for the proposed model 

300 class. Specifically, pointwise confidence bands and bias-corrected confidence 

301 bands are available for the smooth parameter functions; see Nychka (1988), 

302 Ruppert et al. (2003). 

303 The MEM representation of the penalized approach also has the advan- 

304 tage that it can easily accommodate a large variety of effects, at no increase 

305 in the level of complexity of the algorithm itself. In our case, the MEM 

306 representation allows for a unified framework for smoothness selection and 

307 estimation of all model components in (1), including functional random ef- 

308 fects. 

309 4. Identifiability. 

310 4.1. Imposing suitable identifiability constraints. Addditive models for 

311 scalar responses have to ensure identifiability by imposing suitable con- 

312 straints on the functions that make up the addditive predictor (3o + f(x s ). 

313 Most commonly, restrictions like a sum-to-zero constraint Ya=1 fs{ x si) = 

314 for each function f s (x s ) are incorporated in the fit (Wood, 2006a). If this 

315 were omitted, any constant could be added to one function and subtracted 

316 from the others without changing the fit criterion. 

A similar issue arises in the context of our proposed model. For any arbi- 
trary functions g z i(t), Bo(t) a model 

E(Yij(t)) = a(t) + j(zuj,t) + Boi(t) can be rewritten as 

E(Yij(t)) = (a(t)+j zl (t) + B (t)) + {j(z Uj ,t) - 7 zl (t)) + (B 0i (t) - B Q (t)) 

317 to obtain the same fit with a different parameterization. To avoid this, we 

318 impose sum-to-zero constraints for each t so that 

319 n- 1 Y^Lt niB 0i (t) = n- 1 Y,fli E"=i > *) = Vi - 

320 We also center observed trajectories Xf ji j(s) as Xf ^(t) = X f t ij(t) — X f(t) 

321 with mean function Xf(t) = n _1 Ylij Xf,ij(t). If both the sum-to-zero con- 

322 straints for each t are imposed and functional covariates are centered before 

323 entering the model, all effects that vary over the index of the response are 

324 directly interpretable as deviations from the overall mean trajectory a(t). In 
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325 other words, imposing the constraint that the mean value of the effect is zero 

326 at each value of the index of the response function, i.e., Ylil( z i->t) = OVt, 

327 and using detrended functional covariates yields a model that can be in- 

328 terpreted naturally. Standard sum-to-zero constraints implemented in mgcv, 

329 which would correspond to Yli t^( z ^^) = in our representation, yield 

330 neither identifiable models nor effects that are interpretable like this. Im- 

331 plementationwise, we use the method described in Wood (2006a, ch. 1.8.1) 

332 to absorb the sum-to-zero-for-each-t constraints into the design matrices of 

333 all effects varying over t, see section A in the supplement for details and 

334 examples. 

335 4.2. Limits on the identifiability of complex regresssion surfaces for low- 

336 rank functional covariates. For function-on-function-regression terms 

337 J s X(s)(3(s, t)ds, identifiability of /3(s, t) is guaranteed under conditions de- 

338 rived in He et al. (2003), Chiou et al. (2004) and Prchal and Sarda (2007), 

339 which are hard to verify empirically. In practice, an important quantity in 

340 this regard is the effective rank of the covariance operator of X(s), which 

341 is defined as the number of eigenvalues that together account for at least 

342 99.5% of the variation in the covariate. Regression surface estimates can be 

343 unstable if the effective rank of the covariance operator of the functional co- 

344 variates is smaller than the number of marginal basis functions in the tensor 

345 product spline used to represent the coefficient surface (Scheipl and Greven, 

346 2012). 

347 The potential for severe estimation errors is especially large if the kernel 

348 of the functional covariate's covariance operator overlaps the function space 

349 spanned by parameter vectors in the nullspace of the tensor product spline's 

350 roughness penalty. Based on theoretical considerations and simulation re- 

351 suits (c.f. Scheipl and Greven, 2012), we recommend that practicioners check 

352 the effective rank of the observed covariance matrix of functional covariates 

353 and adjust the number of basis functions accordingly. The pffr function 

354 contains such a check and adjusts the number of basis functions if necessary. 

355 5. Empirical evaluation. 

356 5.1. Simulation study. 

357 Simulation Setup. We simulate data for a model Yij(t) = i]ij(t) + £y(t) for 

358 the following four scenarios to investigate the sensitivity of the estimates to 
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359 varying model complexity, noise levels and numbers of subjects and obser- 

360 vations per subject: 

361 1. Functional random intercept and functional random slope: 

362 r]ij(t) = a(t) + B i0 (t) + Bn(t)uij 

363 2. Functional random intercept and one functional covariate: 

364 THj(t) = a{t) + f X Xi ij(8)px(8,t)ds + B i0 (t) 

365 3. Functional random intercept and two functional covariates: 

see rjij(t) = a(t) + / X Mi (y)£i(s,i)ds + / X 2 , ij (s)f3 2 {s,t)ds + B i0 (t) 

367 4. Functional random intercept, one functional covariate, one smooth 

368 scalar covariate effect and one varying coefficient term: 

369 r]ij(t) = a(t) + f Xi )ij (s)^i(s,t)ds + 7i(zi,ij,i) + 8 2 (t)z 2 ,ij + B i0 (t) 

370 Definitions of the various effect functions and descriptions of the data gener- 

371 ating processes used for the covariates can be found in section B of the online 

372 supplement, along with unabridged simulation results and graphical displays 

373 of data and estimated effects for the replications with minimal, maximal and 

374 median error for each scenario. 

375 For each of the four scenarios, we run 10 replications for each combination 

376 of the following settings, yielding 1920 model fits in total: 

377 - number of subjects: M G {10, 100} 

378 - mean number of observations per subject: rtj E {3, 20}. Subject labels i' E 

379 {1, . . . , M} are drawn from a multinomial distribution with probabilities 

380 P(i' = i) oc yi to generate unbalanced designs. 

381 - number of grid points for t: G G {30, 100} 

382 - relative importance of random effects: SNRg E {0.1, 1, 10}, where SNRg is 

383 the ratio of the standard deviation of the linear predictor without random 

384 effects divided by the standard deviation of the random effect functions, 

385 i.e. for SNRg = 1 the contribution of each functional random effect to the 

386 observed trajectories is about equal to the contribution of the non-random 

387 effects, while for SNRb = 10, the contribution of each functional random 

388 effect is about 10 times smaller. 

389 - signal-to-noise ratio: SNR e E {1,5}, where SNR e is the ratio of the stan- 

390 dard deviation of the linear predictor divided by the standard deviation 

391 of the residuals a £ . 

392 Our results show that 10 replications for each combination are sufficient to 

393 derive precise estimates of effects of the setting parameters on estimation 

394 errors and computation times, c.f. Figures 1 and 2. 
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395 Fits are obtained with the default spline bases used in pffrO, i.e., cu- 

396 bic B-spline bases with 20 basis functions and first order difference penalty 

397 for the functional intercept, tensor products of cubic B-spline bases with 

398 five marginal basis functions for the tensor product terms, with first order 

399 difference penalties for the t- and s-directions and second order difference 

400 penalties for the covariate direction (if applicable). The smoothing param- 

401 eters are estimated based on REML as implemented in mgcv. In general, 

402 the models for settings 1 to 4 include between 120 to 1020 spline coeffi- 

403 cients and between 5 to 8 smoothing parameters. Specifically, models have 

404 20 + 10 M coefficients and 5 smoothing parameters for scenario 1, 45 + 5M 

405 coefficients and 5 smoothing parameters for scenario 2, 70 + 5M coefficients 

406 and 7 smoothing parameters for scenario 3, and 75 + 5M coefficients and 8 

407 smoothing parameters for scenario 4. 

408 Simulation Results. We use the relative integrated mean squared error de- 

409 fined as rIMSE(/(f)) = n" 1 £?=i ^JJj{0-^ to evaluate the accuracy 

410 of the estimates. By using the relative errors, we can compare estimation 

411 errors across different scenarios and noise levels regardless of the ranges 

412 of the true ft(t). Note that we are evaluating the estimation accuracy of 

413 the effects on the scale of the response, not on the scale of the coefficient 
4u function itself, i.e., for a functional covariate X(s) we consider the error 

415 f r f s (X(s)P(s,t) - X(s)/3(s,t)) 2 dsdt and not f r J s 0(s,t) - (3(s,t)) 2 dsdt. 

416 A graphical analysis of results (see Appendix B in the supplement) shows 

417 that there are no relevant interaction effects between the setting parameters 

418 M, rii, G, SNRg and SNR £ on the observed errors within scenarios, so we 

419 fit log-linear models with main effects for the setting parameters in each 

420 scenario to observed rIMSE values and proceed to analyze the estimated ef- 

421 fects. Figure 1 shows baseline levels and the estimated multiplicative effects 

422 of the simulation settings on the rIMSE for the estimated quantities. The 

423 effect of increasing the number of gridpoints G for Y(t) from 30 to 100 is 

424 not shown, as it decreased relative errors for all quantities by about half. 

425 Baseline rIMSE values (top left panel) are given for data with SNR e = 1, 

426 M = 10, rii = 3, G = 30, SNRb = 0.1. In this very noisy setting with 

427 small sample size and dominant random effects, covariate effect estimates 

428 are quite bad, with relative errors that exceed one. Since the random effects 

429 are estimated with little error, however, the error for the responses in this dif- 

430 ficult setting is small as well. Increasing SNR e from 1 to 5 (top right panel) 

431 decreases relative errors between fourfold and 16-fold. An increase of the 

432 number of groups from M = 10 to M = 100 (second row, left panel) has no 

433 substantial effect on the overall estimation accuracy of the y(t)-trajectories. 
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Fig 1: Baseline and estimated multiplicative change in rIMSE for the 4 scenarios. The scenarios 
are depicted with different symbols, and the segments accompanying the symbols correspond to the 
estimated effect ± 2 standard errors. Effects other than Bo(t) only occur in a subset of scenarios. 
Horizontal axis on log 2 -scale. 



434 Estimation accuracy of the functional random effects is not improved by this 

435 increase in M either, due to the commensurate increase of the number of 

436 parameters to estimate, while errors for the covariate effects decrease three- 

437 to eightfold. An increase in the average number of observations per group 

438 from rij = 3 to rij = 20 (second row, right panel) has a much larger benefit 

439 on the overall estimation accuracy of the y(t)-trajectories due to the large 

440 improvement of the estimates of the functional random effects. The decrease 

441 is especially strong for scenario 1 which does not include any non-random 

442 effects. The reduction of rIMSE is three- to fourfold in general, and 8-fold 

443 (16- fold) for the functional random intercept (random slope) in scenario 1. 

444 The 7-fold increase of average n$ from 3 to 20 decreases the error between 

445 about three- and sixfold in the same pattern as the tenfold increase in M. 
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446 A reduction of the relative contribution of the random effects to the linear 

447 predictor, i.e. increasing SNRg from 0.1 to 1 [10] (bottom row, left [right] 

448 panel) improves the overall estimation accuracy of Y(t) slightly [two- to 

449 eightfold]. This overall improvement is due to the large reduction of errors 

450 for the covariate effects, which compensates for the observed deterioration 

451 of random effect estimates. While the errors for the former decrease about 

452 four- to 32-fold [four- to 64- fold], the errors for the latter increase 1.5- to 

453 fivefold [15- to 70-fold]. 

454 In summary, our results indicate that estimation accuracy of covariate 

455 effects is affected most strongly by changes in the noise levels SNR e and 

456 especially SNRg, and less strongly by changes in the available number of 

457 observations, i.e. changes in M, rii and 67. Regardless of whether the effect 

458 in question is a simple functional regression coefficient, an index- varying 

459 smooth effect or an effect surface for a functional covariate, the patterns 

460 of relative change in accuracy are identical. The estimation accuracy of the 

461 functional random effects is affected most strongly by the relative importance 

462 of the random effects SNR# as well. While the number of groups has little 

463 effect on the precision of random effect estimates, the effect of the number of 

464 replications per group is strong. In all the settings we considered, important 

465 effects that contribute relevantly to the predictor are estimated with good 

466 to excellent accuracy. Less than 3% of relative integrated errors for Y(t) are 





12h - 




6h - 




3h - 


90 


min - 


45 


min - 


20 


min - 


a 10 


min - 


M 5 


min - 






2 


min - 


1 


min - 




20s - 




10s - 




5s - 




2s - 



enano: 



Computatio n time 

scenario: 2 scenario: 3 



4^ 



scenario: 



: + 



$ 100 



10:3 10:20 100:3100:20 



10:3 10:20 100:3100:20 

M : n. 



10:3 10:20 101)::! 1011:20 10:3 10:20 100:3 100:20 



Fig 2: Computation times for scenarios 1 to 4 (from left to right). Vertical axis on log 10 -scalc. 
Horizontal axis for the various combinations of numbers of subjects M and average number of 
replications per subject 7ij. Results for G = 30 in dark grey and in light grey for G = 100. Timings 
are wall-clock time taken on an 2.2 GHz AMD Opteron 6174. 
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greater than 0.1 - even in the most challenging data situations with few 
noisy observations and small group sizes, our approach is able to reproduce 
the true structure of the data well. 

Computation times. Figure 2 shows computation times on an 2.2 GHz 
AMD Opteron 6174 processor for the different scenarios and sample sizes. 
Especially for models with multiple random effects (scenario 1) computa- 
tion times increase dramatically in M. Smaller models are fit rapidly, and 
even for the largest data sets with N = 2 • 10 5 , computation times are not 
prohibitively long. Speed gains for REML inference on large data sets can 
be achieved by using the pf f r ()-option to use mgcv's bam() routine for es- 
timating additive models on data sets that do not fit into memory, as in 
Section 5.3. Also note that using GCV optimization (Wood, 2004) instead 
of REML-based inference in pf f r can yield up to 10-fold speedups especially 
for large data sets, but tends to be less stable. 

5.2. Predicting Precipitation Profiles from Temperature Curves for the 
Canadian Weather Data. The Canadian weather data consists of tempera- 
ture and precipitation curves, measured as the monthly average over several 
years at 35 Canadian weather stations (see Figure 3, top). The data has been 
used extensively in the functional data analysis literature. As it is available 
as part of the fda-package (Ramsay et al., 2011) for R, we can make our 
analysis fully reproducible (see supplementary material). We will here focus 
on both the functional relationship between temperature and precipitation 
profiles as well as on the spatial nature of the data, clearly visible from the 
locations of the weather stations depicted in Figure 3 (middle left). 

Ramsay and Silverman (2005) propose a concurrent model to predict pre- 
cipitation profiles from temperature curves, where temperature is allowed 
to influence log-precipitation linearly at the same time point t. Within the 
framework of our model, we can investigate more flexible functional regres- 
sion models that allow for lagged temperature effects. Additionally, we can 
take into account the spatial correlation structure between weather stations. 
We consider the model 



Yi(t) = a{t) + 7 ( Ci , d h t) + b 0l + / Xi(s)(3{s, t)ds + en, e it ~ N(0, o*), 



where Yi(t) and Xi{t) denote the log-precipitation and the (centered) tem- 
perature at location i and time t, and Cj and d% denote longitude and latitude 
of location i. As the temperature curves were centered pointwise across sta- 
tions and 7(<3j, di, t) is constrained to sum to zero for each t, a(t) indicates the 



AMMS FOR CORRELATED FUNCTIONS 



19 



495 mean log-precipitation curve for a station with average temperature profile. 

496 The spatio-temporal term j(ci,di,t) yields a smooth cyclic residual curve 

497 for each station, with a spatial covariance structure induced by the bivariate 

498 spline basis for longitude and latitude and its smoothness penalty. Alter- 

499 natively, it can be viewed as a smooth spatial effect that varies cyclically 

500 throughout the year. Small scale local differences in the levels of precipita- 

501 tion not captured by these spatially correlated residuals are modelled with 

502 a scalar random intercept &oi '~ N(0,a%). It is important to note that 



503 the temperature and precipitation profiles considered here are averaged over 

504 several years. Our model thus does not have the problem of 'the future 

505 influencing the past', but relates general weather patterns to one another. 

506 To take into account the cyclic nature of both the response and the predic- 

507 tor curves, we use cyclic basis functions in both s and t direction. The model 

508 can then be fit using the pffrO function in the refund package (Crainiceanu 

509 et al., 2011) as 

pffr(logp ~ s(lat, Ion) + c(s(place, bs = "re")) + 



ff(temp, yind = t, xind = s, splinepars = list(bs = c("cc"))), 
bs.int = list(bs = "cc"), bs.yindex = list(bs = "cc")) 



510 Here, logp and temp are 35 x 12 matrices containing the log-precipitation 

511 and centered temperature profiles, respectively; lat and Ion are vectors con- 

512 taining the latitude and longitude of each location, respectively, and s and t 

513 are vectors containing the months 1 to 12. A smooth overall mean function 

514 a(t) is added by default, s (lat , Ion) yields a smooth spatio-temporal surface 

515 j(ci,di,t). c(s(place, bs="re")) fits a scalar random intercept 6oi f° r the 

516 stations, with c() notation indicating effects constant in t and "re" denoting 

517 a random effect term. A linear functional regression term J Xi(s)(3(s,t)ds 

518 is specified with ff (temp, yind=t, xind=s, splinepars=list(bs=c("cc"))), 

519 with bs=c ( " cc " ) specifying a cyclic cubic regression spline basis in both s and 

520 t direction. The commands bs.int = list(bs = "cc") and 

521 bs.yindex = list(bs="cc") likewise specify a cyclic cubic regression spline 

522 basis for the intercept curve and all other terms that are functions in t, 

523 respectively. 

524 Most of the variation in the data is explained by the effect of temperature, 

525 followed by the spatial effect and the scalar random intercept. The estimated 

526 overall mean function (Figure 3, bottom left) shows a seasonal pattern, with 

527 low precipitation in the spring and high precipitation in the fall. The effect 

528 of temperature on log-precipitation is shown in Figure 3 (bottom right). In 
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529 s direction, a clear seasonal pattern is visible, with higher temperatures in 

530 winter and especially autumn associated with an increase of precipitation, 

531 and higher temperatures in the spring and and summer asssociated with a 

532 decrease of precipitation. Across t, effects are much stronger for the winter 

533 months, and much weaker for the summer months, especially the effect of 

534 winter temperatures on precipitation in the middle of the year; a plausi- 

535 ble result. The spatially correlated smooth residual plus the scalar random 

536 intercept for each weather station is depicted in Figure 3 (middle right). 

537 For a version without overlapping, please see Figure 6 in Appendix B. It 

538 shows some interesting local features, which clearly illustrate that regional 

539 effects cannot capture the spatially varying structure of precipitation curves 

540 adequately. For example, the Arctic station Inuvik in the very north-west 

541 shows a similar error pattern to Continental stations in the north-west (pre- 

542 cipitation higher in winter and lower in the summer and especially in spring) . 

543 These northern Continental stations, in turn, exhibit a pattern which is com- 

544 pletely different than that of the more southern Continental stations (higher 

545 in the summer and lower in the winter). The online supplement contains the 

546 full code for this analysis and explores alternative model formulations for 

547 these data in Section C. 

548 5.3. Application: Modeling fractional anisotropy curves from a DTI study. 

549 The motivating application for this work is a Diffusion Tensor Imaging (DTI) 

550 study comparing the cerebral white matter tracts (WMTs) of multiple scle- 

551 rosis (MS) patients and healthy controls. Fractional anisotropy (FA) is a 

552 function of the three eigenvalues of the estimated diffusion process that is 

553 equal to zero if water diffuses perfectly isotropically (Brownian motion) and 

554 to one if water diffuses anisotropically (perfectly organized and synchronized 

555 movement of all water molecules in one direction) (Heim et al., 2007). 

556 We focus on FA measured at every voxel of several well identified WMTs: 

557 the corpus callosum (CCA) and the right corticospinal tract (CST). Tracts 

558 are registered within and between subjects using biological landmarks iden- 

559 tified by an experienced neuroradiologist. For the purpose of this application 

560 we consider averages of water diffusion measurements over two of the dimen- 

561 sions, i.e. over "slices" through the WMTs along their length, which results 

562 in a functional observation with scalar argument. The study comprises of 

563 106 MS patients and 42 healthy controls, who are observed at one to eight 

564 visits spread over up to 4 years; see Goldsmith et al. (2011), Goldsmith 

565 et al. (2012), Staicu et al. (2011), Ivanescu et al. (2011). Part of these data 

566 is available in the R-package refund. Due to the lack of repeated visits in the 
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Fig 3: Log-precipitation (top left) and temperature (top right) at 35 Canadian weather stations, 
color coded for region. Middle left: Stations' locations. Middle right: Estimated spatially correlated 
smooth residual for each weather station. Bottom left: Estimated overall mean effect. Bottom 
right: Estimated functional effect $(s,t) of temperature in month s on log-precipitation in month 
f, color-coded for sign and pointwise significance (95%): blue if sig. < 0, light blue if < 0, light red 
if > 0, red if sig. > 0. 



22 



F. SCHEIPL ET AL. 



control group, we restrict our analysis to the MS group. Our goal is to bet- 
ter understand the spatial and temporal course of the demyelination process 
in the MS patients. Using only the measurements for the first two visits, 
Ivanescu et al. (2011) investigated the 'spatial' and 'temporal' association 
between the FA along the CCA tract and the FA along the CCA and the 
CTS tracts at the baseline visit for MS subjects. In this paper we investigate 
this aim in more depth by using data from all the consecutive visits reported 
for the MS patients, and by accounting for correlation between repeatedly 
observed functional responses for the same subject. 
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Fig 4: FA profiles along the CCA for three MS patients, one female (red) and two males (blue 
dashed and dash-dotted), at the current visit (left panel). Also depicted are the de-trended and 
smoothed FA profiles along the CCA (middle) and along the CTS (rightmost panel) at their 
previous visit. 



576 Specifically we model the 'temporal' and 'spatial' associations between the 

577 FA along the CCA at the current visit and the FA along both the CCA and 

578 the CTS tracts at the previous visit for consecutive visits, while including 

579 a subject specific random effect. For illustration, Figure 4 shows the FA 

580 profiles at a current visit along the CCA tract for all the 106 MS subjects 

581 investigated (left panel) and the detrended and smoothed FA profiles at the 

582 previous visit along the same tract (middle), and the CTS tract (right). We 

583 consider the model: 

Yij{t) = a{t) + S(t) gi + i{t)zi + v{ Uij )+ 

J X hij (s)px(s, t)ds + J X 2tij (r)fo(r, t)dr + B 0i (t) + e ijt 

584 where Y{j(t) is the FA profile at location t on the CCA tract, Xijj(s) and 

585 X2 ij(r) are the FA profiles at locations s and r along the CCA and CTS 

586 tracts, respectively, observed at the previous visit, i.e., X\^j is the smoothed 

587 and centered previous response Yij-i and e^t ~ N(0, of). Also, m is the age 

588 at the first visit of the ith subject, gi is the gender indicator of the ith subject: 

589 §i = 1 if the ith. subject is male and for female, and Uij is the time lag (in 
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590 years) for the jth visit relative to the first visit of the ith subject. More pre- 

591 cisely, we assume that the mean FA curve consists of two main components. 

592 The first one comprises a gender group-specific function, a(t) for females 

593 and a(t) + 5{t) for males, a varying coefficient function ^(t) for age Z{ to 

594 account for varying effect of age over the tract, and a tract location-invariant 

595 smooth function u{uij) to account for the effect of the elapsed time since the 

596 first visit. The second component captures simultaneously the associations 

597 between the FA on the CCA tract at the current visit and the FA along 

598 the same tract at the previous visit using the bivariate parameter function 

599 Pi(s, t), as well as between the FA at the current visit and the FA along the 
eoo CTS tract at the previous visit through /^(r, t). By taking advantage of the 

601 flexibility of the proposed methods, the variability of the FA profiles at the 

602 current visit is not restricted to be fully specified by the FA profiles at the 

603 previous visit, as is common in functional regression (Ivanescu et al., 2011). 

604 Specifically, subject-specific deviations from the mean FA profile are mod- 

605 eled using functional random intercepts Boi{t). For simplicity of exposition 

606 the random component combines the two types of random effects discussed 

607 in model (1), i.e., Boi(t) comprises the scalar random intercept and the func- 

608 tional random intercept. Because more than half of the MS subjects under 

609 study have only two repeated measurements, we did not include additional 

610 functional random effects such as functional random slopes for the elapsed 

611 time since the first visit. In general, FA is recorded at 93 locations along the 

612 CCA tract and at 55 locations along the CTS tract, so it is reasonable to 

613 assume a dense sample design for both functional predictors. However, while 

614 for most patients the functional measurements are taken at a regular grid 

615 of locations, for several patients measurements are missing at a few loca- 

616 tions along the tract. Furthermore the FA measurements are observed with 

617 noise. Hence, in order to apply the methodology described in this paper, a 

618 pre-processing of the functional measurements that are used as covariates is 

619 needed. Specifically, for each tract in turn, the FA profiles are first detrended 

620 to make the estimated effects easily interpretable (see Section 4.1) and then 

621 smoothed, so that missing values can be imputed. For detrending, the sam- 

622 pie mean curve, which is obtained by averaging the tract-specific FA profiles 

623 over subjects and visits and ignoring the missing observations, is subtracted 

624 from each profile. For smoothing, we employ standard functional principal 

625 component analysis (Di et al., 2009; Yao et al., 2005) to the entire sample 

626 of tract specific FA curves, by using a working independence assumption 

627 between the profiles observed for the same subject. 

628 Then (7) can be fit using the pffrO function in the refund package as 



24 



F. SCHEIPL ET AL. 



pffr(FA ~ c (gender) + gender + c(age) + age + c(s(time_elapsed)) + 
c(s(subj, bs = "re")) + s(subj, bs = "re") + 
f f (FAcca.prev, yind = t, xind = s) + 
f f (FAcst .prev, yind = t, xind = r)) 

629 Here FA is the 254 x 93 matrix of the FA profiles along the CCA at the 

630 current visit; the matrix contains missingness at several locations along the 

631 tract for five subjects. FAcca.prev and FAcst .prev are 254 x 93 and 254 x 55 

632 matrices of the tract specific detrended and smoothed FA profiles along the 

633 CCA and CTS tracts respectively at the previous visit. It is important to 

634 emphasize that, although the functional response involves missingness, no 

635 pre-processing is required; our methodology is applicable, age is a vector con- 

636 taining the age at baseline of the MS subjects (centered and scaled to have 

637 mean-zero and variance equal to one across subjects), and gender and subj 

638 are vectors containing the gender indicator and the subject ID. The smooth 

639 female-group mean function, a(t) is added by default. The specification of 

640 the various model components is analogous to the previous example; for 

641 conciseness, we omit it in this section. Importantly, the overall functional 

642 random intercept Boi(t) in (7) combines a scalar random intercept, pro- 

643 duced by c(s(subj , bs="re")) and a functional random intercept, which is 

644 produced by s(subj, bs="re"). This model explains nearly 83% of the to- 

645 tal variability in the current FA profiles along the CCA tract. Most of the 

646 variation in the fitted model is explained by the functional random intercept 

647 and the effect of the FA-CCA at the previous visit. Figure 7 in Appendix 

648 C shows the various components of the fitted model for 16 randomly drawn 

649 observations. 

650 Figure 5 shows the estimated model components. The estimated mean of 

651 the FA profiles along the CCA tract in women with MS, a(t), is centered 

652 about 0.5 and shows a steep rise at the beginning of the tract, followed by a 

653 sharp decline up to location 20, and continued with a slight increase which 

654 is interrupted by a modest drop at about 70 and then increases rapidly. 

655 The estimated mean function in men with MS, a(t) + 5(t), has almost the 

656 same shape, but with slightly smaller FA values. There are no significant 

657 differences between males and females in FA along the CCA tract. The 

658 estimated varying coefficient function ^/(t) for age at baseline increases in 

659 t and is mostly negative. It indicates that older MS patients tend to have 

660 significantly decreased FA values in the first half of the CCA tract. The 

661 estimated time lag effect, P(-) implies a somewhat cyclic effect of the time 

662 since the first visit: the FA values first decrease in the first 18 months, 
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Gender effect: a{t)(+$(t)) Age effect 7(f) Time effect f>(u) 




n 2<] 40 60 80 
CCA tract location t 



Fig 5: Estimated components of model (7) with ± 2 pointwise standard errors. Top left: mean 
FA in females a(t) (red) and in males a(t) + 8(t) (blue). Top middle: Varying coefficient -y(t) for 
age. Top right: Location-invariant smooth effect 0(u) for lag time relative to the first hospital 
visit. Bottom: Regression functions /3i(s,t) (left) and /3%(r, t) (middle) for the FA profiles along 
the CCA tract and along the CTS tract, respectively, recorded at the previous visit. Color-coded 
for sign and pointwise significance (95%): blue if sig. < 0, light blue if < 0, light red if > 0, red if 
sig. > 0. Bottom right panel: Predicted functional intercepts Bgi(t). 
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663 followed by a considerable increase for the next 18 months, when they begin 

664 to decrease again. The association between the FA profiles along the CCA 

665 at the current visit and the previous visit, $i(s, t), is depicted in the bottom 

666 left panel of Figure 5. Since the association is mostly flat and negative for 

667 the middle region, the results suggest that, for the middle part of the CCA 

668 tract, patients with below population average of FA at the previous visit will 

669 tend to have mostly higher FA at the subsequent visit, whereas patients with 

670 above population average of FA at the previous visit will tend to have mostly 

671 lower FA at the subsequent visit. The positive peak around (s = 0, t = 50) 

672 indicates that below population-average FA at the beginning of the CCA 

673 tract at the previous visit are associated with decreased FA at tract locations 

674 50 or larger at the subsequent visit. Similar interpretations apply for the 

675 peaks around (s = 50, t = 0) and (s = 90, t = 30). The steep positive 

676 peaks around (s = 90, t = 90) and, to a lesser extent, (s = 0, t = 0) show 

677 that there are strong concurrent effects in those regions, i.e. elevated or 

678 decreased FA values in those regions tend to become more extreme over 

679 time. In a similar fashion our results indicate that subjects with FA values 

680 along the CTS tract that are below population average at the beginning of 

681 the tract, but above average for the middle tract tend to have slightly smaller 

682 FA values on the CCA tract at the subsequent visit. Overall, the effect of 

683 FA-CST is much smaller than the effect of FA-CCA at the previous visit. 

684 As a caveat, note that significance here is pointwise significance, so that 

685 some of the observed features may not hold up in a simultaneous assessment 

686 of significance. Simulation results also indicate that estimation accuracy of 

687 fixed effects may be low in settings in which the random effects dominate the 

688 predictor such as this one. Finally, the predicted intercept curves i?oi( - ) 

689 displayed in the bottom panel. The large variation in the predicted functional 

690 intercepts reveals large inter-subject variability. 

691 In conclusion, modeling the FA profiles along the CCA tract using our 

692 flexible modeling framework shows that most of the variability in the data is 

693 captured by subject-specific random effects. Accounting for the FA profiles 

694 at the previous visit can provide new insights into the spatial and temporal 

695 development of the disease. In particular our results point to mostly weak 

696 associations between the current FA values along the CTS tract and the 

697 subsequent FA values at the CCA tract. 

698 6. Discussion and Outlook. We propose a comprehensive framework 

699 for flexible additive regression models for correlated functional responses, 

700 allowing for multiple functional random effects with, for example, spatial 
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701 (Section 5.2), temporal, spatio-temporal or longitudinal (Section 5.3) cor- 

702 relation structures. Dependence structures can be modeled either implicitly 

703 by specifying smooth temporal, spatial or tempo-spatial effects or explicitly 

704 by including functional random effects with correlation structures that are 

705 constant over the index of the response and are given by the precision ma- 

706 trices of Gaussian (Markov) random fields. Inference in this framework can 

707 be performed by standard GAM software, allowing us to take advantage of 

708 established robust and flexible estimation algorithms. The approach is im- 

709 plemented as fully documented open-source software in the pf frO-function 

710 in the refund package (Crainiceanu et al., 2011) for R. 

711 Simulation experiments show that the proposed method recovers relevant 

712 effects reliably and handles small group sizes and/or low numbers of repli- 

713 cations well. Data sets of considerable size can be fit in acceptable time 
7u with the proposed method. Two application examples demonstrate that the 

715 model class we propose makes it possible to fit flexible model specifications 

716 that do justice to complex data situations and yet still yield interpretable 

717 results that help understand the underlying processes. 

718 This work opens up a number of interesting avenues for further research. 

719 One direction is to extend the type of relationship between the functional 

720 covariates and the response. Specifically, the current model setup implies 

721 that all values of a functional predictor influence the response trajectory 

722 across its whole range. However, the proposed model can easily be modified 

723 to accommodate fixed, limited integration ranges, e.g. to allow for cumu- 

724 lative effects f se u hl t _/j 2 i X(s)(3(s, t)ds over intervals of fixed lengths for 

725 temporal covariates and responses that are observed on the same domain, 

726 where responses cannot be influenced by future covariate values as suggested 

727 by Malfait and Ramsay (2003). In the limit, this also includes the simple 

728 concurrent model with terms X(t)/3(t). An implementation of this feature 

729 will be included in a future version of pffr(). 

730 A second direction regards improving the estimation of the covariance struc- 

731 ture of the residuals. Since the inference algorithms we use at present do not 

732 exploit the extreme sparsity of the design matrices for smooth observation- 

733 specific residual terms, estimating such terms dramatically increases compu- 

734 tation time and memory requirements for large data sets. On the other hand, 

735 simply assuming i. i. d. errors en will often be unrealistic since some degree 

736 of auto-correlation and heteroscedasticity over the index of the functional 

737 response is usually encountered in practice. An iterative procedure similar 

738 to the approach in Reiss et al. (2010), where observed residuals from an ini- 

739 tial model estimated under a working independence assumption are used to 
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740 estimate a working covariance structure and the model is then re-estimated 

741 based on de-correlated data may offer an efficient alternative approach that 

742 we plan to investigate in the future. 

743 In a third direction, we are currently developing diagnostic measures to iden- 

744 tify potential problems with low-rank functional covariates (c.f. Section 4) 

745 as well as practical model-building strategies regarding the estimation of 

746 corresponding regression surfaces. 

747 Further on, implementation of a dedicated toolbox for REML-based infer- 

748 ence tailored to function-on-function regression based on the computation- 

749 ally efficient array regression approach of Currie et al. (2006) may facilitate 

750 inference in large scale problems and help to generalize the proposed meth- 

751 ods for multivariate functions and image regression. 
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Fig 6: Spatially correlated smooth residual plus scalar random intercept for each weather station. 
Stations are roughly ordered from north-west to south-east within regions. Solid lines for the 
spatially correlated random effect, dots for the observed residuals. Color coding is as in Figure 3. 




Fig 7: Detrended data and estimated model components for the DTI application for 16 randomly 
selected observations: Grey lines: Observed, detrended data Yij(t) — »(i)(solid line), detrended 
fit Yij(t) — d(t)(dotted line); Blue lines: predicted functional random intercepts Boi(t); Black 
lines: estimated contribution of functional covariates f X\ij (s)$i (s, t)ds + f X2 } ij(r)02{i",t)d.r; 
Red lines: estimated contribution of scalar covariates (male gender, age, time lag since first visit) 
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